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A comprehensive set of motion equations for a multi-actuated flight vehicle is presented. The dynamics 
are derived from a vector approach that generalizes the classical linear perturbation equations for flexible 
launch vehicles into a coupled three-dimensional model. The effects of nozzle and aerosurface inertial coupling, 
sloshing propellant, and elasticity are incorporated without restrictions on the position, orientation, or number 
of model elements. The present formulation is well suited to matrix implementation for large-scale linear 
stability and sensitivity analysis and is also shown to be extensible to nonlinear time-domain simulation through 
the application of a special form of Lagrange’s equations in quasi-coordinates. The model is validated through 
frequency-domain response comparison with a high-fidelity planar implementation. 


I. Introduction 

Modeling the behavior of a flexible booster in atmospheric flight is a non-trivial exercise in dynamics that has been 
treated with various levels of fidelity over the last fifty years by practitioners of the art. For the purposes of basic 
sizing and performance analysis, it is often sufficient to assume the vehicle is a rigid body or even a point mass and 
propagate the trajectory considering as few as two degrees of freedom. Traditionally, physical vehicle and trajectory 
design considerations have helped to ensure that the simulated behavior was closely replicated in the actual system; 
for example, great expense is incurred in eliminating nonlinearity from the physical system so as to ensure that the 
flight control stability problem is tractable using classical linear techniques. 

Despite the intentional design simplifications, simulation and analysis of this class of flight vehicle requires a 
carefully constructed and well-understood dynamics model that considers the significant coupling effects between 
the rigid body dynamics, structural modes, sloshing propellant, aerodynamics, servodynamics, and the flight control 
system. Control- structure interaction is often the chief design challenge, especially in boosters with large slenderness 
ratio and low structural mass. 

In an era of emerging unconventional approaches to the design of launch vehicles, missiles, and hypersonic aircraft, 
traditional (namely, planar) dynamics models 1,2 may no longer be sufficient for modeling the vehicle dynamics for the 
purposes of stability analysis and control design. Traditionally, by virtue of a set of body axes generally aligned with 
the principle moments of inertia, it could be assumed that the control system axes were fundamentally decoupled. 
Clearly, this is not possible in the more general case, and more importantly, the structural response may have a strong 
multi-axis character. Finally, multi- actuation is not well-treated by traditional models, which may assume a set of 
identically performing, non-interacting actuators and nozzle masses that are lumped together in order to compute 
inertial coupling effects. * 1 

The intent of the present formulation is to overcome these and other restrictions in fidelity and scalability while 
retaining the general character of the linear perturbation models from which this new approach originates. Indeed, 
the present formulation is derived using almost the same technique as in the planar case, albeit much more rigorously 
owing to the matrix-vector nature of the mechanics. Ultimately, the result is a compact recipe for modeling flexible 
aerospace systems with multiple attached effectors, including provisions for sloshing propellant. With parameters fixed 
in time, the equations of motion are readily implemented in linear form so as to generate a linear state-space model 
for control design. If the nonlinear terms are included, the formulation can be integrated numerically as a high-fidelity 
mechanization of the full vehicle dynamics. Thus, the application of the model is twofold: rather than propagate 
a nonlinear trajectory simulation and load the resultant parameters into an analytic stability model (or numerically 
linearize a nonlinear simulation), the nonlinear equations of motion can be directly linearized at certain operating 
points for analysis. 

The present formulation is known as an integrated formulation, which differs from a multi-body formulation. In 
a multi-body formulation, the equations of motion for each body are written with respect to an inertial frame, and 
the necessary interbody forces are solved as a set of adjoint or constraint equations that relate the motion of each 
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body to its adjacent bodies. 3 The constraints can be eliminated analytically to give a minimal realization, or they 
can be eliminated algorithmically. 4 Such formulations are appropriate when the structure of the system undergoes 
large kinematic displacements, but are intractably complex if one seeks to generate a compact, explicit description 
of the system motion equations. Notably, some semi-rigid multibody codes are based on an open tree topology, 5 as 
is common in modeling spacecraft appendages attached to a rigid bus. In contrast, the present system is also of tree 
structure, but consists of a flexible base body to which rigid actuated appendages are attached. 

In the case of the flexible rocket vehicle, it is convenient that the attached effectors (nozzles, aerosurfaces) undergo 
only small rotations from their nominal positions, and the effector mass is small with respect to the total system 
mass. As a result, the dependency of the integrated system moment of inertia on the motion can be neglected, and the 
primary motion coordinates are selected to be those of the center of mass of the entire system. The primary coordinates 
are often called quasi-coordinates , 6,7 since they are not inertial coordinates but instead expressed with respect to the 
instantaneous axes of a rotating frame. 

The dynamics of the integrated vehicle are derived in vector form using kinematic assumptions that allow linearity 
to be retained in the expressions of kinetic and potential energies. The approach could be considered a vector extension 
of the classic linear perturbation equations of motion for a rocket, 1,8 and the dominant nonlinearities (gyroscopic 
effects) in the resultant motion can be included via application of a special form of Lagrange’s equations. 6,7,9 However, 
one intent of the present approach is the generation of a three-axis coupled linear stability analysis model, which does 
not require the inclusion of the nonlinear terms. 

The paper is organized as follows: Section 2 defines the physical quantities and notation used in the derivation and 
formally states assumptions related to the underlying mechanics. Section 3 presents the derivation of the equations of 
motion and the output relations in complete form. Section 4 provides an overview of the approach to mechanization of 
the motion equations in a scalable numerical code; simulation results are presented. Section 5 concludes the discussion 
and summarizes the results. 


II. Definitions and Assumptions 

The flexible rocket system is modeled as a set of interacting subsystems. The scope of the present model is limited 
to generalized forces due to thrust and aerodynamics; all other forces arise internally. The system inputs will be 
assumed to be a set of effector torques r ei that rotate each effector about an appropriate pivot. In a rocket system, 
this torque is typically provided by a hydraulic actuator, 10 and the resultant generalized force on the system is angle- 
dependent, resulting from a deflection of a thrust vector or an aerodynamic normal force on the effector itself. 

System outputs consist of the vehicle truth states, augmented by a set of output equations that model the states as 
detected by a suite of sensors at particular locations fixed in the body frame. Sensed vibration effects are of paramount 
importance in feedback control of flexible structures; the sensor output models ensure that the physics are captured 
with adequate fidelity. 

The relevant elements of the vehicle configuration are shown in Figure 1 . The integrated vehicle is partitioned into 
a base body B , to which a frame of the same designation is rigidly attached, a set of r effectors, and a collection of n 
linear oscillators. The B frame rotates and translates with the base body. The remainder of the dynamics are described 
as relative motion coordinates, whereby motion of elements other than the base body degrees of freedom are written 
with respect to B. Note that B is arbitrarily located in space, although we will later designate the origin of B to coincide 
with the integrated vehicle center-of-mass. 

The base body itself is assumed to be a rigid body of mass m v with an inertia tensor J v about B , located with position 
r v with respect to B. 

Each nozzle or aerosurface is considered to be itself a rigid body; in each effector is embedded one of r gimbal 
frames G/, the motion of which with respect to B is described by the angular rate co Gr Each effector has mass m ei and 
inertia tensor J ei about Gi and is located at displacement r ei with respect to each G* . Each gimbal Gi is displaced from 
the origin of B by the vector r Gi . 

The sloshing fluid in propellant tanks is modeled using the usual mechanical analog 11 where, for small pertur- 
bations, the liquid involved in sloshing undergoes linear motion and is partitioned from the liquid assumed rigidly 
affixed to the vehicle structure. The sloshing component of the liquid is represented by a spring-mass damper system 
vibrating in a plane normal to the thrust axis; the sloshing natural frequency and damping are a function of the axial 
acceleration and the container geometry. Each j of n slosh masses is a point mass of magnitude m s j described by its 
unperturbed location r sj with respect to B\ its relative position vector p sj describes its local translation. The inertia of 
the sloshing fluid is assumed invariant with respect to liquid motion and is given by J S j. 
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Figure 1. Flexible vehicle with effector and sloshing mass 


Finally, the flexibility of the system is assumed to be in mass -normalized diagonal form; i.e., the orthogonal 
modes of a finite element model. The finite element model with p degrees of freedom is described by a vector 
generalized displacement rj e R p , generalized velocity fj e R p , and a set of associated partial eigenvectors £ 

M? xp corresponding to the translational and rotational components of the eigenvectors, respectively, at each discrete 
gridpoint. The p degrees of freedom contain only the vehicle elastic modes; it is assumed that the rigid-body modes 
(zero-eigenvalued modes) are both orthogonal to the elastic modes and have been removed from the finite element 
model via an appropriate partitioning procedure. 

A summary of the variables describing the system is shown in the following table. 


Quantity 

Description 

B 

Body frame, fixed with respect to rigid base body 

p Gi 

Location of gimbal frame of each i of r effectors with respect to B 

m v , J v , r v 

Mass, moment of inertia, and location of vehicle CM with respect to B , 
not including effectors and sloshing propellant 

Wlei ? J eU p ei 

Mass, moment of inertia, and location of of i th effector with respect to Gi 

Wlsj-> Jsj, p sji Psj 

Mass, inertia, position, and displacement of each j of n slosh masses 


Elastic generalized displacement and generalized velocity 


Table 1. Dynamical variables 


Fundamentally, we assume that certain tensors in B are unchanged with respect to angular displacements. For 
example, the value of an inertia tensor of the i th effector about frame G, coordinated in frame B, varies with the 
direction cosines that relate the basis vectors of B and G. The variation is assumed negligible for small rotations. If 
the two frames are nominally aligned, 

Jei = VC BG V^[C GB ] * J G . (1) 

In assuming small rotations, the transformation effecting a rotation between two frames A and B is sequence-independent 
and can be partitioned into a linear form such that 

[C BA ] = [/ + £*], (2) 

where I is the identity and s x is the skew-symmetric form of the small rotator. We will also use the skew- symmetric 
form to represent the vector cross-product. 6 
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In the following section, we derive the equations of motion for the integrated system. 


III. Dynamics 

The approach to the derivation of the dynamics yields both linear and partial nonlinear forms. First, the kinetic and 
potential energies are described in terms of mixed velocities and displacements. Thereafter, the generalized momenta 
are derived for the composite system. Finally, a state vector is chosen and a special form of Lagrange’s equations in 
quasi-coordinates is used to derive the equations of motion. The linearized equations of motion follow as a special 
case. 

A. Kinetic Energy 

The kinetic energy of the integrated system is expressed in mixed coordinates. The kinetic energy of the base 
body is expressed in terms of the translational and angular velocities of the B frame, augmented by the energy of the 
effectors and slosh masses resulting from their local velocity with respect to B. We will assume for convenience and 
without loss of generality that (1) the finite element associated with each effector has its origin coincident with each G> 
frame and (2) the vehicle and sloshing masses are point masses, to which each has associated a single physical elastic 
displacement but an arbitrary number of elastic coordinates. The latter assumption avoids unnecessary bookkeeping of 
integrals or sums associated with elasticity in either continuum or discrete form, respectively, and makes no difference 
in the final formulation. 

The kinetic energy T is given by the sum of the contributions of the empty vehicle and fixed propellant, sloshing 
propellant, and effectors, such that 


T = T v + T sj + Tei- (3) 

7-1 i= 1 

The velocity of the vehicle base body mass in the inertial frame is 

v v = v B + 0,7) + cu%r v (4) 

where e R 3xp . The velocity consists of the reference frame velocity, the velocity due to elastic motion, 1 " and the 
velocity due to rotation of the B frame. Likewise, the velocity of a slosh mass in the rotating frame is 

v sj = v B + <3> sj ri + oj x r sj + p sj , (5) 

and the effector velocity is 

v ei = V B + Og/7) + ^ B (r Gi + r ei ) + cjQfei + Q¥ G ifj) x r ei , (6) 

where we have accounted for the additional angular velocity of each frame with respect to B. 

It follows then that the elements of (3) can be expanded as 


2T V 

= m v (u H + <£> v f] + co B r v ) 2 


2 Tsj 

= m s j (v B + O s jf] + aj B r s j + p s jf 

(7) 

2 T ei 

= Wl e i {v B + ®GiH ^ B (VGi *'ei) + + ^GiH) ^eij 



where the operator (.) 2 in the vector context is taken to be the appropriate transpose-inner-product operation. 

Tt is helpful to clarify that the orthogonality of the finite element modes is with respect to motion of the integrated vehicle, not the base body. 
Therefore, the expression (4) describes the sum of elastic motion of all mass excluding the effectors and slosh masses. 
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B. Generalized Momenta 


1. Rigid Linear Momentum 

From the kinetic energies, the generalized momenta of each subsystem can be calculated. The generalized mo- 
menta with respect to body frame translation are determined to be 

W VB T v = m v v B + m v <b v fi + m v aj^r v 

= m s jv B + m S j<t> S jfj + m S jO>^r S j + m sj p s j (8) 

V VB T e i = m ei v B + m ei Q> Gi f] + >n el co y H (r Gi + r ei ) + m ei co Gi r ei + m ei ^ Gi i]) x r ei 

The terms in (8) are grouped such that 

V Vb T v = m v v B - m v r x co B + m v <S> v f] 

V„ B r s; - = m S jVB-m S jr*u>B + m S jp S j+m sj ® sj ri (9) 

V DB T ei = m ei v B - m ei {r Gi + r ei ) x u B - m ei r x Lj Gi + - m ei r x A! Gi fp 

note that if m v + Z m S j + Z m ei - (total mass) and that if we now assume that B is at the center of mass of the 
integrated vehicle, 

n r 

X m A r A + Z mei(l ' Gl + r «') + m v r » = °> (1°) 

7=1 *'= 1 

which, as expected, decouples the translational momentum from the angular momentum since the velocity ^ is that 
of the center of mass. 

We note that momentum is conserved in vibration by virtue of orthogonality, so 

n r 

Ws/h,// + - m ei r x >¥ Gi i]) = 0, (11) 

7=1 ( =1 

and the total linear momentum of the integrated body is 

n r 

V Vb T = m T v B + ^ m s jp s - ^ m ei r x co Gi . (12) 

7=1 1=1 

The last term in (12) is the change in linear momentum due to angular velocity of an attached effector and is the 
lateral effect of “tail-wags-dog” inertial coupling. 

2. Rigid Angular Momentum 

The generalized momentum of the body frame with respect to the angular velocity ojb is found to be 
V Ws T v = m D r x v B + m v r x O v f] + m v r x 0 J B r v 

Vco B T sj = m S jr* j v B + m sj r*<& sj Ti + m S jr*u B r S j+m s r x p S j (13) 

V Wg T ei = m ei (r Gi + r ei ) x v B + m ei (r Gi + r ei ) x O Gi f] 

+m ei (r Gi + r ei ) x (o B (r G i + r ei ) 

+m ei (r Gi + r ei ) x 0 J Gi r ei + m ei (r Gi + r ei ) x n , Gi f]) x r ei . 

Thereafter, we group terms and again note that angular momentum is conserved in normal vibration and that the 
origin of B is the center of mass. We also employ the identity 

L m i r f 0J B r i = Jl °B- ( 14 ) 

N 

The resultant simplifications yield 

— d v (x>B 

Vaj B T S j = JsjOJB + msjr^jpsj (15) 

^ cjs^ei ~ Jei^B + WleiiXGi + r ei) ^Gi^ei’ 
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The effector energy can be expanded such that 


V WB T ei = JeiOJB - ^yPiei^ojGi + m ei r y a) Gi r ei . 

The last term, using (14), is the effector inertia about the gimbal pivot point; 

^ ei ~ Jei^B + (7 ei ~ ? Gi^ed e \)^Gi' 


(16) 


(17) 


Then, letting the integrated system inertia be J T = J v + 2 J S j + Z Jeu the total angular momentum with respect to 

B is 

n r 

Vw b T = J t ^b + 2 m sj r*jp sj + J ei - r Gi m ei r y )a> Gi . (18) 

1=1 <=i 

As in the translational case, the last term in (18) is change in angular momentum of the rigid body due to the 
angular velocity of an attached effector, or the angular effect of “tail-wags-dog.” 


3. Momenta Due to Slosh 

The generalized momenta with respect to the slosh coordinates are found to be 

= 0 

\m sj v B + m sj <& sj ri + m sj co*r sj + m sj p sj , j = a 
\ 0 , jtcr 

= 0 . 

The total momentum of each sloshing mass is therefore 

Vp s T = m sj p sj + m sj v B + m sj O sj fj - m sj r x sj co B . 




Tscr 


V T 

v psj 1 el 


4. Momenta Due to Effectors 

The generalized momentum of each effector due to its angular velocity is given by 


VovTv = 0 

= 0 


V T 

v a>Gi 1 ecr 


m ei r x v B + m ei r y ^> Gi f] + m ei r y oj y (r Gl + r ei ) 
+m ei r* j u> Gi r e i + m ei r^^¥ Gi f}) x r e i, 

0, 


l = (T 

j 4 (T. 


Reordering the terms in (21), it follows that 


Va> Gi Tei = Jei^Gi + m ei r x v B + m ei r x co B (r Gi + r ei ) + m ei r x ^ Gi f] + m ei r^ G ihYr ei \ 
noting that m e ;r*uj*(rG; + r ei ) = (J ei - meir^r^cos, the angular momentum becomes 

V Wg .T = J ei (x) Gi + m ei r x v B + ( Jei ~ m ei r y r Gi )CL> B + (m ei r*® Gi + Jei'Vcdfj. 


5. Momenta Due to Bending 

The momentum due to bending can be determined as 

V /; 7’„ = <b v m v v B + <b v m v <& v f] + <b v m v u) B r v 

V'/ T SJ ^ sjWlsjVB "1" ®sjWlsj®sjf] "1" ^ sjWlsjGJ B r s j + Q>sjm s jPsj 

Vi,T ei = f ei m e iV B + t ei m e i^ G if] + t ei m ei a) B (r G i + r ei ) 
+^ei m ei a> G i r ei " 1 " ^ ei^ei(^ GiP) ^ei 
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(19) 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


( 24 ) 


where 


r ei = 0 Gi - r%¥ Gi . (25) 

This term contains the rotary inertia of the attached nozzle. Assuming again that modal coordinates are orthogonal 
to rigid body coordinates and that the body frame is fixed at the center of mass, 


Vf,T v = 

® v m v <& v f] 


^f/Tsj = 

jfn.y + o, jltl s jp s j 

(26) 

VfjTei = 

teimei^Gif] + t eitrieiLJg^ei + f G ,-7?) X r ei . 


Note that the generalized mass can be defined as 


1 u = 0„OT„0„ + 

n r 

^ ^ &sjWl sj^sj + ^ ^ {pei^ei^Gi ~ f ei^eF ' e i^ Gij ? 

(27) 


7=1 i= 1 


which includes the rotating nozzle masses. Finally, expand the remaining term in (26) as 

i^ei^Qi^e i ~ ^ Gi^ei ^ei^Gi ^ei ~ ^Gi^eF e i^Gi ~ i^Gi^ei ~ ^Gi^ei^ e ij ^Gi ( 28 ) 

and the bending generalized momentum becomes 

n r 

~ Ed + ^ ^ ^ 5 jMsjPs j + ^ v i^GiJei ~ ^Gi^ei r e ij tO Gi • (29) 

7=1 i= 1 

The last term in (29) is the contribution to the bending momentum due to the nozzle angular velocity. 

C. Potential Energy and Damping 

The system potential accounts for the strain energy due to bending, the potential due to deflection of each slosh 
mass, and the potential due to the torsional stiffness of each effector rotated around its hinge point. Likewise, the 
energy dissipation is a function of the bending generalized velocity, the slosh perturbation velocity, and the angular 
rate of each effector. The expressions for bending are derived from the well-known second-order linear structural 
dynamic model recast in a modal form via a similarity transformation x = ^>rj such that 

OMOf] + D b rj + OTCO// = O/ (30) 

with OAIO = yu, <D7C® = K b , and some assumed damping matrix D b . While generally ji - 7, it is not strictly 
necessary that D b ,K b be diagonal or have any particular properties. It is often necessary to admit asymmetric damping 
and stiffness that account for aeroelastic coupling effects. 12 As can be seen from the developments in Sections (1) and 
(2), it is only required that the elastic response is orthogonal to the rigid-body degrees of freedom. 

It follows that the potential and dissipation functions can be written as 

n r 

2V = f]K b rj + ^PsjKsjPsj + ^ 0 G iK ei 6 Gi (31) 

7=1 i= 1 

n r 

2 D = fjD b fj + ^ ] psjDsjPsj + ^ v F) G iD e iCO G i. 

7=1 i= 1 

The assumption has been made that f o Gi dt ~ 0 Gi ; kinematics of the nozzle have been ignored on account of small 
rotations. 

D. Generalized Forces 

7. Generalized Forces Due to Effector Motion 

The primary exciter of bending, lateral, and angular motion is the generalized force due to thrust or control surface 
moments. Suppose there exists an effector force in the body frame modeled as a unit vector of undeflected force u ei 
and its associated magnitude R a . The force in the body frame of the i th effector is 

fei = \ C BG ]R ei Uei (32) 
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where [C BG ] is the kinematic transformation that expresses the rotating unit vector u ei in the body frame. If the 
rotations are small with respect to the nominal cant angle, the transformation [C BG ] can be partitioned; 

fei = [l + &Gi + CM*] Re i Uei = Rei^ei ~ ReiKfioi ~ ReiU*'¥ Gi*l ■ (33) 

The resultant moment is then 


T ei ~ r Gifd ~ r Gi^ei u ei r Gi^ei u e i^Gi r Gi^ei u e i^ Gi6- 


(34) 


Clearly, there is a static component and a component that is dependent on the generalized coordinates. The ex- 
pressions (33) and (34) are suitable for modeling a thrusting, gimbaled rocket nozzle. Note that we do not complicate 
the dynamics be constraining the rotation 6a to occur about a particular pair of axes; since the nominal cant angle 
u e i can be arbitrary in the body frame, a rotation in the body frame is also, in general, a three-axis rotation. For a 
rocket nozzle, servo torques are applied in the appropriate local pitch-yaw sense (“rock” and “tilt”) via a constant 
input transformation that maps the local gimbal coordinate system into the body frame. The nozzle stiffness K ei and 
damping D ei are adjusted to approximate the constrained nozzle along the local roll axis in the gimbal frame. 

If the effector type is a simple aerosurface, it can also be modeled as a linear angle-dependent moment effector, 
where R e i is an effectiveness matrix that varies as a function of flight condition. 


2. Generalized Forces Due to Aerodynamics 

It is straightforward to add a complex nonlinear aerodynamics model using the typical table-lookup method to 
apply forces to the rigid body degrees of freedom. For stability analysis, it is often necessary to consider perturbations 
with respect to a trim condition. In the case of an accelerating rocket, the reference frame is taken to be the accelerating 
frame tangent to a gravity turn trajectory. The angle of attack of the vehicle is a function of the vehicle pointing angle 
6b and lateral velocity v T with respect to the trajectory. For small angles, the transformation from body to trajectory - 
relative acceleration is given by 

vt = vb + 6^as (35) 

where as is the vector of instantaneous thrust acceleration expressed in the body frame. The equations of motion are 
then transformed to trajectory frame relative motion by inserting the inverse transformation. Thereafter, the linear 
aerodynamic force and moment on the rigid body are approximated by a set of aerodynamic stiffness and damping 
matrices, where 


faero ~ ~ ( + C tt Vr + Q r (Z>g) 

^aero ~ ~ i^rr^B + CftVj 1 + C rr (i)B ) • (36) 

The relevant terms embedded therein are a function of dynamic pressure and the 6-DoF linear stability derivatives. 


E. Boltzmann-Hamel Equations 

The equations of motion of the integrated system can be determined from the generalized momenta and a modified 
form of Lagrange’s equations. Note that the energy has already been expressed in the most convenient form for 
analysis. Let us write Lagrange’s equations in the classical form for an unconstrained system; 


d [dL\_ dL 

dt \ dqi ) dqi 


(37) 


where L = T - V, q t are each generalized coordinate and Q t is the generalized force associated with that coordinate. 
For this special class of system, we note that d T / dqi = 0; it is sufficient to rewrite (37) in vector form as 


d 

dt 


(v 4 r) + v q v = q. 


( 38 ) 


Without loss of generality, let us choose a state vector of true velocities 


v 

Cu 

O77 
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where v are holonomic, or integrable velocities, and u are non-holonomic or quasi- velocities. The transformation C 
is the orthonormal kinematic transformation that coordinates the non-holonomic velocities u in the inertial frame, and 
the transformation <D is the orthonormal linear transformation that maps elastic generalized velocities 77 into physical 
velocities. 

The direct vector form of Lagrange’s equation applies to the holonomic velocities, and observing that there are no 
position-dependent potential terms, 


4 (V„D = Q u . (39) 

at 

In the case of the quasi- velocities u , we substitute the transformed coordinates into (38) and find that 

cicV u T) + CV u T + V f V = Q u (40) 

at 

where £ = Cu. Recall that C relates an inertial to a rotating frame and is always orthonormal, so CC = /. The time 
derivative C = Ccj* and equation (40) can be written as 

J t (y u T) + 0)gV u T + CV f V = CQ U . (41) 

Finally, we consider the elastic velocities’ contribution to the kinetic and potential energies. The elastic physical 
velocities are orthogonalized via some constant linear transformation. The elastic velocities are indeed relative motion 
coordinates with respect to the rotating frame, but they are coordinated with respect to an inertial basis. The elastic 
coordinates are therefore holonomic coordinates. Physically, the structural displacements must coincide with the 
undeformed rigid body, which follows from assumption that the rigid body modes are to be excluded from the elastic 
solution. 

Since Lagrange’s equations for the elastic coordinates undergo only a constant transformation, we have that 

®Jt ( v,?r ) + ^ ” v = Qn (42) 

and, with O® = /, the flexibility equations have the expected classical form 

J t (v,^) + (43) 

The expressions (39), (41), and (43) are a special vector form of the Boltzmann-Hamel equations, or Lagrange’s 
equations in quasi-coordinates . 9 They are much simplified compared to the more general forms of the Boltzmann- 
Hamel equations. 
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F. Equations of Motion 

The time-invariant equations of motion in final form follow directly from the generalized momenta and the appli- 
cation of the Boltzmann-Hamel equations. Under the previously stated kinematic and dynamic modeling assumptions, 
the complete trajectory-relative motion equations are given below. 



In the full nonlinear form given above, the motion equations are easily integrated into a trajectory simulation via 
partitioning into the requisite generalized mass matrix E. In time-invariant form, the linear components are readily 
mechanized into a first-order descriptor form by writing the linear dynamics as 

Ex = Ax + Bu (49) 

y = Cx + Du 


with the state vector 

[ Ob 0GI 0q2 ... 0 Gr Psi p s 2 . . . psn V 1 1 

X= r _ _ _ _ — — T J T nr • (50) 

[ Vb Cl>b Cl>G\ MG! • • • <^>Gr Psl Ps2 • • • Psn ?7 J 

Linear models of the flight control system and the servoactuators are block-diagonalized with the plant dynamics 
model and the multiple-input, multiple-output (MIMO) frequency response and eigenvalues can be computed numer- 
ically using commonly available computational tools. 
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G. Output Models 

The system outputs capture the effects of sensor mounting on an arbitrary point on the structure. We assume 
that the sensor output axes are nominally aligned with the body axes. Sensor types to be considered in the present 
application are rate gyros, integrating gyros (attitude angle sensors), and accelerometers. 

Each attitude angle and rate sensor detects the superposition of the body attitude with the corruption due to elastic 
motion, given by 

Om = Os+Vyrnri (51) 

~ ^ ymJ] 

where e M? xp is the subset of column eigenvectors associated with rotation at the m th sensor. 

The accelerometer sensor includes the linear effects of local rotation and the nonlinear effects due to non-colocation 
with the center of mass. If r ym is the vector from the center of mass to the sensor in the body frame, the sensed 
acceleration with respect to the reference frame is 

Yym ~ [G ] (vb + p ym + ^^B P y^ ^ B P y fn + ^ b^B P • 

Note that due to flexibility, the transformation [ C YB ] rotates the sensed body acceleration through small angles; 
also, the sensor velocity and acceleration with respect to the center of mass are a function of the elastic generalized 
coordinates. We neglect the small velocity of the sensor relative to the center of mass due to burning propellant. It 
follows that 

Yym — ^^yrnjf) J {^B + . (52) 

In most applications, we can safely ignore the higher-order flexibility terms and model the accelerometer as 

yym ~ VB + ( j J B r y>n + ^ B^ B P y fn ~ p B^ymJl + 2£Ug0^ m ?) + (53) 

IV. Simulation 

The present formulation has been incorporated into a general missile autopilot stability analysis and design tool, 
FRACTAL (Frequency Response Analysis and Comparison Tool Assuming Linearity). FRACTAL is a scalable im- 
plementation of the equations of motion that allows the trajectory to be simulated at a various flight conditions under 
the assumptions of fixed (linear, time-invariant) or time-varying parameters so as to assess the impact of parametric 
uncertainty in a given vehicle configuration. While the current mechanization now allows full nonlinear trajectory 
simulation in the software, FRACTAL’S primary analysis task is the generation of time and frequency domain results 
to characterize the stability characteristics of the integrated system using linear control system stability metrics. 

Owing to the easily implemented matrix-vector form of the equations of motion, numerical assembly of linear sys- 
tems representing the full vehicle dynamics with autopilot and actuators in-the-loop is exceptionally rapid on modern 
computing hardware. The software is equipped to automatically construct the equations of motion to arbitrary scale up 
to the limits of memory and computational capacity; for example, a typical FRACTAL linear model of a five-engine 
heavy-lift launcher concept may have as many as 260 states. The efficient form of the equations of motion makes 
the model particularly well suited to numerical optimization of the flight control system design parameters, which 
employs extensive finite differencing implemented in a parallel computing environment. 

Validation 

Various simulation test cases have been constructed to compare the response of the present formulation with well- 
established and validated dynamics models. A time-domain test case example is shown below in Figures 2 and 3, 
where the FRACTAL model dynamics are a trajectory-relative perturbation model. The results have been compared 
to TREETOPS, 13 a multi-body dynamics package based on Kane’s equations. 3 In the simulation results shown, an 
example system consisting of a rigid vehicle with an attached nozzle is assembled together with a pair of linear 
servoactuators that provide torques about the nozzle in the pitch and yaw axes. The vehicle gimbal angle command is 
excited with a small-amplitude sinusoidal wavelet to verify the coupled response of the thrusting engine, the actuators, 
and the vehicle to which they are attached. With small angles and small body rates, gyroscopic coupling is minimal 
and the models agree exceptionally well even considering the omitted nonlinearities. 
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TWD test case 1 : body attitude 




Figure 2. Time domain simulation, pitch attitude response 




time (sec) 

Figure 3. Time domain simulation, pitch rate response 


A second simulation test case is shown in Figure 4. In this simulation, a special parameterization of the FRACTAL 
multi-axis simulation model is compared with the high-fidelity planar equations of motion. 2 In both simulations, the 
nozzle, actuator dynamics, flexibility, slosh, aerodynamics, and sensor output models are included in computing the 
frequency response from the gimbal angle command to a rate gyro affixed to the vehicle structure. The simulations 
show excellent agreement, confirming that the planar dynamics are indeed a special case of the present formulation. 



-4 -3 -2 -1 0 1 2 3 

10 10 10 10 10 10 10 10 

Frequency (Hz) 

Figure 4. Frequency domain simulation, pitch rate response 


V. Discussion 

The work detailed in this paper has illustrated a technique that generalizes the planar linear perturbation dynamics 
models that have been the staple of launch vehicle control design efforts for over five decades. The resultant formula- 
tion is based on a careful balance of assumptions that intends to retain the maximum possible fidelity while eliminating 
the complexity that lends itself to only a strictly algorithmic implementation. By writing the system energies in a dom- 
inantly linear form, the salient features of the true physics are retained and the higher-order and relatively insignificant 
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kinematic effects have been omitted. Most importantly, the formulation leads to an explicit and intuitive set of motion 
equations that are easily coded in both linear and nonlinear forms. The model has direct applications to Monte Carlo 
simulation, stability analysis, and numerical optimization, and is extensible to a wide variety of vehicle configurations, 
including launch vehicles, missiles, and thrust- vector-controlled aircraft. 
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Motivation 


• Future launch vehicles to be considered for development are not planar, 
single- or lumped-nozzle launch vehicles (LVs) 

• Ascent flight control design and dynamic modeling for asymmetric systems is 
a challenging design task (for example, Shuttle) 

• Vehicle may be inertially, elastically, and aerodynamically asymmetric, with 
multiple nozzles/effectors of varying types 

• No easy-to-use, integrated simulation environment exists for modeling these 
systems that can be simultaneously used for flight control design and stability 
analysis as well as high-fidelity nonlinear dynamic simulation 

• Existing multiple-degree of freedom multiple-body dynamic simulation tools 
are not always the best choice 

• Classical perturbation equations are high fidelity, but only for planar and 
inertially symmetric systems 
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Flight Dynamics Features of Large LVs 


• Thrust vector and aerosurface control 

• Canted, asymmetrically located nozzles 

• High flexibility due to low structural mass 

• Large potential for control-structure and fluid-structure 
interaction 

• Large, slow actuator systems with inertial and servoelastic 
coupling 

• Large sloshing masses with lightly damped response 

• Large and unpredictable bias forces and torques 

• Inertial cross-coupling of control axes 

• Time-varying mass properties 

• Time-varying environment 
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A New Formulation 


• We propose a new flight mechanics formulation based on a vector extension 
of the classical perturbation equations for launch vehicles 

• Nozzle and aerosurface inertial coupling, sloshing propellant, and elasticity 
are generalized in the model 

• Derived from a rigorous energy-based approach that leverages well-founded 
assumptions to minimize the formulation’s complexity 

• Results in a closed-form set of motion equations 

- A nonlinear simulation can be formed from Lagrange’s equations in 

quasi-coordinates (the Boltzmann-Hamel equations) to propagate the trajectory 

The linear equations can be directly used to form a state model of the system 
dynamics for control design 
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Modeling Assumptions (I): Overview 


• The vehicle is divided into a base body (B), gimbaled effectors (G), and 
sloshing propellant (s) 

• Position and orientation of the vehicle elements are arbitrary 


m 

SJ 
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Modeling Assumptions (II) : Integrated Vehicle 


• External generalized forces are due to thrust and aerodynamics 

• Model inputs are a set of effector torques that rotate nozzles/aerosurfaces 

• Model outputs are truth states and a set of sensor outputs at arbitrary 
locations on the body frame 

• All elements are coordinated in the body frame and undergo small 
rotations/displacements (nonlinear kinematics are negligible) 

Launch vehicle-specific assumption; not well-suited to spacecraft with flexible 
appendages or deployable masts 

Change in integrated MOI due to effector motion is small 

• Vehicle mass properties vary slowly enough that higher-order effects due to 
element time derivatives are ignored 


MOI = moment of inertia 
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Modeling Assumptions (III): Effectors 


• Each effector is a rigid body attached to a 
local gimbal frame with some angular velocity 

• The effector has a 3x3 MOI, stiffness, and 
damping matrix, as well as at least one plane 
of actuation and a nominal thrust direction 
along some unit vector 

• Each effector is attached to a single FEM 
gridpoint via a 3-axis hinge; the thrust 
structure near this point is quasi-rigid 

• The mass/inertia and some of the local 
backup structure flexibility of the nozzle 
“pendulum mode” is included in the FEM; for 
small angles, the total deflection (force 
follower effect) is found via superposition 



MOI = moment of inertia 
FEM = finite element model 
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Modeling Assumptions (IV): Slosh 


• Sloshing propellant is modeled using the standard 
mechanical spring-mass-damper analogy 

• The sloshing fluid is a point mass vibrating in a 
plane normal to the mean axial acceleration 

• The stiffness and damping are nonlinear functions 
of the axial acceleration and container geometry 

• The liquid not involved in sloshing is retained in the 
rigid-body mass and moments of inertia 

• The tank is quasi-rigid and is coupled to the elastic 
model by equal force distribution of centerline grids 
along the tank axis 



Energy | Environment | National Security | Health | Critical Infrastructure 



From Science to Solution 


©Science Applications International Corporation. All rights reserved. SAIC and the SAIC logo are registered trademarks of Science Applications International Corporation in the U.S. and/or other countries. 


© co 



Modeling Assumptions (V): Elasticity 


• The elasticity of the integrated system is modeled as a set of linear orthogonal 
modes of a FEM; rigid-body modes are removed via partitioning 

• If the FEM is provided at discrete times, the non-orthogonality (with respect to 
the rigid body) of mass integrals at other mass configurations is considered 
negligible 

• The FEM stiffness and damping matrices need not be symmetric, positive 
definite, or have any particular properties 

• Aeroelastic effects are integrated directly via modification of the stiffness at 
damping matrices as a function of flight condition 


FEM = finite element model 
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Modeling Approach (I): Energy 


• The total kinetic energy of the vehicle is given by the sum of the energies of 
its elements, 

n r 

T = T v + ^ T s j + ^ T ei 

;=] i= 1 

v v = v B + 0> 0 T] + co B r v 
v s j — ug + O s jtj + co B r s j + p s j 
L-ei — t-B + + r e j ) + (x)Q-r e i + an) *ei 


• The contribution to kinetic energies of the elements is due to angular velocity 
of the vehicle, flexibility, and local motion 

• Higher-order flexibility and nonlinear terms are negligible and are omitted 
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Energy Comparison 


Planar versus vector formulation 


Planar Lagrangian 

1 . \2 1 

2 


r=^m^ z +/ v cp + (// v n) 2 +^m45 z +5 s +/ 5 (p + (/f i J7) 2 +|/w n (<5 z -(/ 5 + / fJ )(p+^ & /7-/„(^+^ g ^)r 


Vehicle 


Slosh 


Nozzle 


Vector Lagrangian 


SJ 

OT . 

JLI ei 


2 T v = m v (v B + <\\f] + OJ X B r v f | Vehicle 


2 T sj — fn s j ^ sp] s j P s jj | Slosh~ 


= m ei (v B + O Gi fj + a)g(r G i + r ei ) + (x)%r ei + C¥ Gi fj) x r ei J 


Nozzle 
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Modeling Approach (II): Generalized Momenta 


The generalized momenta follow from the kinetic energy: 

n r 

V 0s r = rn T V B + 2_J m sjPs - X m * T ri a) Gi 

/= l /= 1 

« r 

m sj r *jPsj + 2_j( J ei - r Gi m ei r ri) a> Gi 

j= I *'=> 

V Pl .r = /Hjjps; + m sj v B + m S ji[> S jf] - m sj r*jO) B 


+ m ei r*v B + (Jei - + (/««,,/■£ <P Gi + Jei'i'ciM 


n 


^ jjT — JJt] + ^ ^ ^ ^ f/ 0(5 i 7?X e //’ ei ) 

>1 i'=l 


The time-derivative of the momenta can be readily assembled into the 
generalized mass matrix 

Each term is readily associated with a physical phenomena 
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Modeling Approach (III): Thrust Forces 


• Generalized forces arise due to deflection of a 
thrust vector: 



[ C BG ] R ei u ei 


/ + *£,+ Wcm)' 


Roiit 




— R e jU e i R e iti e -0(ji R e iii e ^Y Qii] 



T ei ~ r Gifei 



• Note that the elastic component can be folded 
back into the stiffness matrix. 

• 2-DoF nozzle motion (rock/tilt) can be reduced 
via the use of a suitable input transformation 
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Modeling Approach (IV): Nonlinear Motion 
Equations 


• The nonlinear motion equations can be readily obtained from the linear 
generalized momenta expressions in the rotating frame and a special form of 
Lagrange’s equations in quasi-coordinates, or the Boltzmann-Hamel 
equations. 

• For this system, 


d / dL 



dt \dcfi / 




+ V q V = Q 


• Choose the state vector 



V 


Holonomic velocity 

q = \ 

\ Cu 


Quasi-velocity 


[ <$77 j 


Elastic velocity 
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Construction of Nonlinear Motion Equations 


• For the quasi-velocities, we use f = Cu and find that 

(V U D + cv^r + V f y = q u 

at 

• C is an orthonormal kinematic transformation, so CC = / and C = 

^ (v„r) + + cv f y = eg. 

• For the flexibility equations, we have 



d 

dt 


since cl) is constant. 
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Motion Equations Example: Angular Motion 


• The full nonlinear equations can be realized via simple substitution. 

• For example, the angular motion equation, neglecting aerodynamics and 
other external forces is 



+ ^ ] Rei ( r Gi Uei r Gi 11 ei * r Gi 11 ei ^ G i ^ ) 

1 

Thrust effect 
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Motion Equations Example: Slosh Motion 


• The application of the special form of Lagrange’s equations yields the correct 
expressions for the acceleration of the slosh mass in the rotating frame 



• More straightforward, compared to a Newton-Euler approach 
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Modeling Approach (V): Sensor Outputs 


• Sensor outputs are formed from rigid body “truth” plus flex corruption 

- Angle and gyro outputs contain the angle component due to flex 


0 m — Og + ^ ymO 

&Bm ~ ^ tjmV 

Accelerometer expression can be derived as 


y ~ + ?ym + ^.COgFym + OJgFy m + OJgOJgFyf^ 

Excluding high-order terms, a sufficient accelerometer model is 

yym ~ &B + ^ ^ ~ ^ ymO + 2cOg0)y m i] + 
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Model Implementation and Verification Activity 


• The present model has been integrated into the generalized dynamic 
simulation, autopilot design, and stability analysis tool FRACTAL (Frequency 
Response Analysis and Comparison Tool Assuming Linearity). 

• FRACTAL began as a planar perturbation tool but now implements the full- 
featured nonlinear dynamics with an analytic linearization capability. 

• Analytic linearizations can be used in parallelizable optimization algorithms for 
flight control gain and bending filter design routines. 

• High-order (>250 state) linear systems can be automatically generated with 
all actuator dynamics and autopilot dynamics included. 

• Verification test cases were constructed against the previously, rigorously 
exercised classical planar dynamics. 

• For special parameterizations, the three-axis model can be made to simulate 
a planar vehicle as a special case. 

• Time-domain verification against Kane’s equations was also performed. 
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Verification Results (I) 






Rigid body + actuator, nozzle inertial coupling, aerodynamics only 
For planar, inertially diagonal linearized system, duplicate response is 

eXpeCted Plant ra ^ e response from pitch gimbal command 
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Verification Results (II) 


• All dynamics (rigid, flex, slosh, slosh coupling, nozzle-flex coupling) 


Rant rate response from pitch gimbal command 
From: noz 1r ock To: Out(1) 



Frequency (Hz) 
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Summary and Forward Work 


• The classical planar equations of motion, based on verified assumptions for 
launch vehicle dynamics, have been extended into a fully coupled, nonlinear, 
multiple degree-of-freedom set of motion equations. 

• The compact form of the final motion equations is easily coded. 

• The motion equations can be analytically linearized without resorting to 
numerical perturbation methods. 

• The code is suitable for any flexible flight vehicle with thrust-vector or 
aerosurface control and sloshing propellant, including launch vehicles, 
aircraft, and missiles. 

• Forward work will examine higher-fidelity time-domain simulation test cases to 
quantify the loss of fidelity due to kinematic assumptions. 

• The present model is being used as a foundation for additional research into 

high-efficiency control allocation algorithms and stability theory for flight 
vehicle adaptive control. m s 
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